Precursors, aftershocks, criticality and self-organized criticality 
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We present a simple model of earthquakes on a pre-existing 
hierarchical fault network. The system self-organizes on long 
time scales in a stationary state with a power law Gutenberg- 
Richter distribution of earthquake sizes. The largest fault car- 
ries irregular great earthquakes preceded by precursors devel- 
oping over long time scales and followed by aftershocks obey- 
ing an Omori's law. The cumulative energy released by pre- 
cursors follows a time-to-failure power law with log-periodic 
structures, qualifying a large event as an effective dynami- 
cal (depinning) critical point. Down the hierarchy, smaller 
earthquakes exhibit the same phenomenology, albeit with in- 
creasing irregularities. 

PACS numbers: 72.10.Fk, 73.40.Hm, 75.10.Fk. 



Seismologists model dynamical rupture propagation 
with complex friction laws and barriers and attempt 
to ascribe the earthquake complexity to nonUnear pro- 
cesses and/or heterogeneities. From a more global point 
of view, it has been suggested that earthquakes are some- 
what similar to critical points and can be addressed 
using tools of the renormalization group. A very broad, 
but quite ill-defined, perspective is also available with the 
concept of self-organized criticality ||,^. These various 
points of view model different properties at different time 
scales; it is hard to see how they relate to each other, and 
whether they are part of a unique meaningful "theory of 
earthquakes" . A closely related puzzle is whether criti- 
cality and self organized criticality are compatible. 

The present work attempts to unify a significant frac- 
tion of this earthquake phenomenology, and to answer 
this puzzle: we define a simple model that exhibits the 
self-critical organization of the crust at large time scales, 
the critical nature of large earthquakes and the short- 
time rupture dynamic properties. 

We start with the following ingredients, (i) The faults 
are organized in a hierarchical geometrical structure [^,^ . 
We do not address the problem of the construction of 
the fault patterns themselves which involves much larger 
times scale (10^~^ years) compared to the time scales we 
describe (10°~^ years), (ii) The tectonic plate is driven at 
a slow average uniform rate and we take into account its 
heterogeneities and the existence of relaxation processes 



by allowing for fluctuations in the local rate of loading. 
{in) When a threshold is reached, a redistribution occurs 
on adjacent faults, with amplitude controlled by the size 
of the faults. 



FIG. 1. The fractal cell structure. 

Our model is an hybridization of the sandpile model |^ 
and of the fractal automaton . As in , the cell sizes 
are arranged as a discrete fractal lattice. Each cell can be 
viewed as representing the region which is elastically un- 
loaded when a fault fails. The fractal distribution of cell 
sizes then represents a fractal distribution of fault lengths 
The number of fractal generations does not appear 
to be a crucial parameter. Most simulations were carried 
out with 8 generations, but some runs with 12 gener- 
ations did not exhibit major differences. We load the 
system by dropping particles at regular intervals (which 
we use as a clock) onto the grid at random sites. The ad- 
dition of a particle is analogous to energy loading. The 
probability that a particle is added to a particular cell 
is proportional to the area A of this cell. Although the 
exact cell (fault domain) to receive the next increment in 
stress is random, the entire grid is loaded uniformly at 
a uniform rate over the long-term. This represents the 
long-term uniform strain at the boundaries between mov- 
ing tectonic plates. The short term random heterogeneity 
in loading represents heterogeneity in crustal structure 
or in upper mantle flow and the associated relaxation 
processes. Each cell becomes unstable when it contains 
n X 4j4 particles, n being a parameter. It then breaks and 
redistributes AA of its particles to its immediate neigh- 
bors. The number of particles redistributed to an adja- 
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cent cell is proportional to the linear dimension of the 
cell. The (n — l) x A A particles that are not redistributed 
are considered as lost, like the particles which are re- 
distributed outside the grid on the plate border. Such 
energy losses we call "cooling". Since the particles rep- 
resent energy, the model assumes that a fault fails when 
the stored energy reaches a critical threshold. The key 
difference with ||] is that the energy must reach the crit- 
ical level over the entire area of a cell before it is allowed 
to break. Due to the fractal structure, cells of widely 
different sizes are thus coupled together, mimicking the 
multi-scale interactions between faults. 

The clock is defined by the particle drops. Cascades 
are triggered by the addition of a single additional par- 
ticle, i.e. they occur instantaneously (a delay can also 
be introduced in the aftershock sequences, see below). 
At variance with the rules of and in accord with the 
standard sandpile model [D , the size of an earthquake is 
determined by the size of the cascade, and is proportional 
to the total number of particles cascading. We thus iden- 
tify the cascade as complexity in the mainshock, e.g. the 
linking of fault strands and segments or even the linking 
of adjacent faults such as occured in the 1994 Landers 
earthquake We define precursors as those regional 
cascades which precede a mainshock and the aftershocks 
as the cascades which follow it. We then use the time 
scale defined by the particle drops to explore the tem- 
poral structure of both the foreshock and aftershock se- 
quences. 

The cooling, i.e. the disappearence of particles dur- 
ing an event, represents a loss of stored elastic energy 
due to the earthquake. Consider for instance an elastic 
medium under constant applied shear strain, which sud- 
denly undergoes a rupture in the form of a dislocation or 
a crack. The stress field is redistributed with enhance- 
ment at the crack tip as well as screening at some other 
places, while both the total shear stress and total elastic 
energy decrease. The amount of loss is a function of the 
nature of the rupture and of geometry. For instance, in 
the Griffith problem, a crack of length 2c is introduced 
into a rod of length L and section ~ L^. Under con- 
stant strain, the relative energy loss ^ incurred by the 

introduction of the crack is ^/{l + ^) @ with a 
depending on the geometry and of the order of unity. If 
we take c between L/4 and L/2, we get ^ anywhere 
between about 0.2 — 0.6. This loss corresponds to fric- 
tion on the fault plane, the creation of surface energy in 
the extensive crushing which occurs in the fault zone , 
and in the radiation of elastic waves |10| whose energy 
is ultimately lost as heat. Comparing with ^ = in 
our model, we see that n ~ 2 is not unreasonable. We 
take this value for most of our simulations below. These 
large uncertainties in cooling are actually not critical - 
in fact, provided n is not too close to 1, the results are 
largely independent of n. The choice n = 1 of |0] leading 



to internal conservation (except at the boundaries) does 
not seem relevant for earthquakes. 

As in 10], we identify an event cascading through the 
largest cell on the grid as a main shock for our system, 
i.e., the largest regional event. After a transient depend- 
ing on the initial conditions, the system self-organizes in a 
stationary state with a power law Gutenberg- Richter dis- 
tribution of earthquake sizes P{E)dE - E-^'^-^^^-^'^dE. 
This can be called a self-organized critical state, mea- 
sured by a statistics encompassing many times the largest 
time intervals between the largest earthquakes. We gen- 
erated about 100 main shocks on an 8*'* order fractal 
structure by adding a total of about 4.6 ■ 10^ particles to 
the system. The average number of time units (particle 
drops) between main shocks is T = 4.6 • 10^, with fluc- 
tuations of this interval of the order of 10%. This quasi- 
periodicity occurs because our main shocks are charac- 
teristic earthquakes ||ll[ which completely control the en- 
ergy redistribution at large scales. Their existence is 
not in contradiction with the self-organized critical state: 
they simply arise because of finite size effects. 

For the purpose of comparison with real earthquakes, 
we choose our units of time so that T = 100 years, i.e. 
roughly 10* time units correspond to two years. For most 
main shocks, precursory built-up of activity is clearly vis- 
ible, and lasts 1 — 2-10^ time units, or about 20 — 40 years. 
A decaying activity posterior to the main shocks is also 
visible with a lifetime of about 0.1 — .2 • 10^ time units, 
i.e. about 2 — 4 years. Quite often, the interval between 
two successive main shocks will not be as quiet, but is 
interrupted by the breaking of one of the second largest 
cells. The time interval between such smaller shocks is 
much less regular than the time interval between main 
shocks. 

We start by discussing aftershocks. The common wis- 
dom holds that aftershocks involve a time delay between 
the application of the stress and the subsequent rupture. 
This delay presumably involves an intermediate relax- 
ation time, whose effect could resemble that of diffusion 
or visco-elastic processes. In the present model, the spa- 
tial heterogeneity of the loading rate already refiects the 
existence of delay mechanisms. Inspired by the analysis 
in terms of critical phenomena (see below), we then plot 
the cumulative number of events as a function of time |l^ ] 
starting from an initial date of about 10* time units pos- 
terior to tc and going backwards in time. If Omori's Law 
^IT^ oc holds, this gives n{t) ^ {t - t^f-P (the 

theoretical divergence at long times for p < 1 is trun- 
cated by the existence of background seismicity). Fits 
performed for 15 of our model aftershock sequences gave 
a distribution of the exponent p centered around p = .9, 
with small fluctuations, p G [0.85,1.05], in good agree- 
ment with the exponent measured for earthquakes. These 
results are independent of the value of the dissipation pa- 
rameter n, provided n is not too close to 1. As n — > 1, 
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in particular when n = 1, as in Q], the exponent p is 
often near 1, although it has large fluctuations, and can 
be found as low as p — 0.7. Simulations on other frac- 
tals indicate that the exponent p is also near unity for 
different geometry. 

The model can be improved by the addition of a delay 
mechanism with a characteristic time t for the avalanche 
associated with the main shock (this has no qualitative 
effect on the statistics of events, nor on the analysis of 
precursors). This had the drawback of introducing one 
additional parameter, but allows a more natural analysis 
of the data. Fortunately, results were found totally in- 
sensitive to the value of r over a wide range (r < 10'*) 
and in perfect agreement with Omori's law again. 
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FIG. 2. Cumulative number of afttershocks as a 
function of time (in units of 10^ particle drops). The 
full line is a power law with 1 — p — .12. 



We now turn to precursors. While virtually all large 
shallow earthquakes have easily recognizable aftershock 
sequences, the same can not be said for foreshocks. In 
fact, if foreshocks are defined to occur with the same 
time and space clustering as aftershocks, then most large 
events do not have a recognizable foreshock sequence 
It is only when the time scale is extended to tens of years 
and the space to hundreds of kilometers that precursor 
sequences can be recognized p^p^. We thus use the 
term "precursor" to distinguish the two definitions. We 
analyze data in the same spirit as in ||l6|,^ : wc plot the 
cumulative Benioff strain e(t) (square root of the energy 
release in an event) as a function of time, starting back 
in time within a range 1 — 3 ■ 10^ (20 — 60 years) prior 
to the main shock. We use a lower cut-off such as to 
exclude the crowd of small background events and test for 
various cut-off values. The conclusion is, in most cases, 
a reasonably clear evidence for a power law behavior, 
decorated by log periodic oscillations ||lq,p^ 



e{t) = A- B{t^ - <)'" [1 + Ccoscjln(ic - t)] 



(1) 



This example is quite typical. Much better fits are some- 
times obtained. Bad fits were rare but sometimes oc- 
cured. The improvment of the when using (|l|) com- 
pared with a pure power law fit (C = 0) was always 
greater than two. The power law in (|l|) is the signature 
of a critical dynamical behavior, reminiscent of depinning 
transition. 
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The fit by this equation of the cumulative foreshock time 
series for the biggest main shock is shown in Figure 3. 



FIG. 3. Plot of Benioff strain before the biggest 
shock. The dotted line is power law behaviour, full line 
is power law decorated with log-periodic oscillations. 
The time is in units of 10^ particle drops. 



The log-periodic correction in (^ reflects the discrete 
scale invariance (DSI) of the hierarchical fault network 
on which the events occur [ pTf (we observed similar oscil- 
lations decorating the Gutenberg-Richter law). Intrigu- 
ingly, the amplitude C of these oscillations is much bigger 
than for equilibrium statistical mechanics models [ p^ : 
the threshold dynamics in the earthquake model seems 
much more sensitive to DSI. We discuss below how the 
value of CO is related to our fractal structure. 

Values of m and lu fluctuate more than the exponent p 
for aftershocks. In 80% of the cases however, m was found 
in the interval m £ [0.2,0.6], in good agreement with 
experimental data |]l6| , ^ . ui was found in the interval 
u! g [6, 12] corresponding to 1.7 < A < 2.8 where u; = 
27r/ log A. 

These results are not sensitive to the dissipation pa- 
rameter n. Only the time scale is modified. In the ex- 
treme case of n = 1 (without cooling), the seismic activ- 
ity is much more random, and the rate of events looks 
much more constant between main shocks than it does 
with cooling: only about 10"* time units (2 years) before 
the main shock does the cumulative Benioff strain de- 
velop a power law behavior on the approach to the main 
shock, with properties that are then similar to the case 

Beside the Benioff strain we studied the correla- 
tion length ^ defined as the maximum spreading distance 
of the cascades. To get more stable numerical results, 
we calculated its integrated value, which does exhibit a 
power law singularity described by 
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C oc (tc - 1)-"^ , t<tc, i^ie [0.6, 0.8] 
^cx(<-^J-'^^ t>t„ (2) 

i/i and V2 have smaller fluctuations than the exponent m 
of the Benioff strain. Observe that they are different on 
both sides of the critical point, a situation which is known 
to be possible in disordered systems in particular. The 
divergence of ^ in confirms the critical point picture. 
Moreover, the exponent provides a relation between DSI 
in time and DSI in space. From our fractal geometry, the 
latter is characterized by a; — > 2a;. Substituting ^ — > 2^ 
in (|) implies \t-tc\ {X)''^\t-tc\ with A G [2.4-3.2], 
i.e u) € [5.4,7.2], in reasonable agreement with what we 
directly observed. 

We now use the time-to-failure law (||) to try to fore- 
cast the main shock by fitting it to the "experimental" 
data up to a cutoff time prior to the main shock, as pro- 
posed in Overall, wc find that 95% of the main 
shocks can be predicted with an uncertainty less than a 
year, four years in advance. This must be compared with 
the typical fluctuations of about 10 years of the time in- 
tervals between two main consecutive shocks. There are 
however cases where the predictability is much higher, 
and also extreme cases where the precursory activity is 
essentially nonexistent. 

We also studied events occuring on the second largest 
cells. The analysis is slightly more complicated because 
two such events can be close in time but well-separated in 
space. We therefore restricted our attention to well iso- 
lated cases. The analysis of precursors and aftershocks 
gives similar results as for the foregoing events on largest 
cells. The only difference is that the time scales involved 
are shorter (roughly by a factor of 10), and the fluctua- 
tions in exponents somewhat bigger. This is presumably 
due to the fact that the relative size of the fluctuations 
of the "energy" field with respect to that of these earth- 
quakes is larger for smaller events due to the infiuence of 
the earthquakes at the upper levels. 

Our results thus indicate that reasonable intermediate- 
time earthquake prediction may be achievable, as pro- 
posed in |l^,|lj]. From the simulations, there is clear 
evidence that the predictabilty depends on the "temper- 
ature" of the system - the larger the loss of energy after 
a main shock, the better is the prediction of the next 
one. This is particularly clear for the first main shocks 
obtained when initiating the empty system, for which 
predictabilty is very high. The astonishing accuracy ob- 
served in [ p^ for the Loma Prieta example might thus 
be due to the existence, in that case, of a very "cool" 
seismic system. 

In addition, we have demonstrated for the first time the 
possible coexistence of self-organized criticality and criti- 
cality. Up to now, they were considered as dual (mutually 
exclusive) modes of behavior: critical depinning occurs 
when the applied force reaches a critical value beyond 
which the system moves globally, while self-organized 



criticality needs a slow driving velocity and describes 
the jerky steady-state of the system. The critical na- 
ture of our large cascades emerges from the interplay 
between the long-range stress-stress correlations of the 
self-organized critical state and the hierarchical geomet- 
rical structure: a given level of the hierarchical rupture 
is like a critical point to all the lower levels, albeit with a 
finite size. The finite size effects are thus intrinsic to the 
process. 
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